function [auxOut,errorMes] = DSGEmodel_yieldCurve_ss(params)

%% Section 1: 
% The size of the model
ny    = 49;     %Number of control variables
nx    = 6;      %Number of state variables
ne    = 5;      %Number of shocks
nx1   = 1;      %Number of endogenous state variables
mx    = 0;      %Number of endogenous state variables (must come first in x)
myx   = 1;      %Number of lagged control variables appearing as states. 
                %These lagged control variables must come first in y AND must
                %appear in x after the endogenous states. 

% Setting the errorMes (1 for a problem, else 0)
errorMes = 0;

%Unfold parameters
U0       = params.U0;
U0d      = params.U0d;
BETTA    = params.BETTA;
B        = params.B;
CHI      = params.CHI;
CHI0     = params.CHI0;
THETA    = params.THETA;
DELTA    = params.DELTA;
ALFA     = params.ALFA;
PHI      = params.PHI;
ZI       = params.ZI;
KAPAw    = params.KAPAw;
NU       = params.NU;
ETA      = params.ETA;
RHOR     = params.RHOR;
PHIpai   = params.PHIpai;
PHIpai_1 = params.PHIpai_1;
PHIy     = params.PHIy;
PHIy_1   = params.PHIy_1;
PHIc     = params.PHIc;
PHIc_1   = params.PHIc_1;
PHIl     = params.PHIl;
RHOz     = params.RHOz;
RHOn     = params.RHOn;
RHOd     = params.RHOd;
RHOp     = params.RHOp;
RHOa     = params.RHOa;

MUZss    = params.MUZss;
PAIss    = params.PAIss;
lss      = params.lss;
KoY      = params.KoY;


% The eta-matrix
eta = zeros(nx,ne);
eta(2,1) = params.SIGMAmuz;
eta(3,2) = params.SIGMAd;
eta(4,3) = params.SIGMAn;
eta(5,4) = params.SIGMAp;
eta(6,5) = params.SIGMAa;


% The higher order moments
% M.M2 is the expected value of kron(eps,eps). 
% M.M3 is the expected value of kron(eps,eps,eps) and so on. 
% Here I assume that the shock is standard normal.
% M.M2=1; M.M3=0; M.M4=3; M.M5=0; 
% If the shocks are independent standard normal you can use the command:
momEps =gaussian_moments(ne);
%% Section 2: Solving for the steady state
% The real stoch discount factor
Mreal      = BETTA*MUZss^(-CHI*(1-CHI0)-CHI0);

% The nominal interest rate
Rss        = PAIss/Mreal;

% % The one-period nominal bond price
% P1ss        = 1/Rss;

% % The one-period real bond price
% P1Realss    = Mreal;

% Value of capital
Kss        = ((4*KoY)^(1/(1-THETA)))*lss;  

% The output level
OUTPUTss   = Kss^THETA*lss^(1-THETA);

% The consumption level
Css        = (1-ZI/2*(PAIss/PAIss^NU-1)^2)*OUTPUTss - DELTA*Kss;

% Marginal costs
MCss       = 1/(ETA*OUTPUTss)*(ZI*(PAIss/PAIss^NU-1)*OUTPUTss*PAIss/PAIss^NU ...
               -(1-ETA)*OUTPUTss - ZI*Mreal*(PAIss/PAIss^NU-1)*PAIss/PAIss^NU*OUTPUTss*MUZss);

% The real wage
Wss        = MCss*(1-THETA)*Kss^THETA*lss^(-THETA);

% The parameter on disutility
PHIzero    = ((Css-B*Css*MUZss^-1)/Css^CHI0)^-CHI*Wss*(1-lss)^(1/PHI)/Css^CHI0;

% The Value function (minus)
VFss       = -(U0+U0d+1/(1-CHI)*((Css-B*Css*MUZss^-1)/Css^CHI0)^(1-CHI) + PHIzero*(1-lss)^(1-1/PHI)/(1-1/PHI))/(1-BETTA*MUZss^((1-CHI)*(1-CHI0)));
AA         = VFss;
EVFss      = (VFss*MUZss^((1-CHI)*(1-CHI0))/AA)^(1-ALFA);

% Constraint on the effective discount factor
if BETTA*MUZss^((1-CHI)*(1-CHI0)) >= 1
    errorMes = 1;
end

% The yield curve
i             = 1:40;
P             = (1/Rss).^i;
nameP = cell(1,40);
for i=1:40
    nameP{i} = ['p',num2str(i),'_t'];
end

%% Section 3: Reporting the output for the perturbation approximation
% For level approximation: k_cu = K;
% For log transformation: k_cu = log(K);
% For logistic transformation: k_cu = -log(1/Vss-1)

%The level of the states
%                    c_ba1  ,muz_cu ,d_cu  ,n_cu   ,paistar_cu  ,a_cu
auxOut.Xss        = [Css    MUZss    1      1         1        1 ]';                      
auxOut.transformX = [1      1        1      1         1        1 ];     
auxOut.labelx     = [{'$c_{t-1}$'},{'$\mu _{z,t}$'},{'$d_{t}$'},...
                     {'$n_{t}$'},{'$\pi^*_{t}$'},{'$a_{t}$'}];



%The level of the controls
%                   c_cu   r_cu   pai_cu  w_cu    l_cu   output_cu  mc_cu  evf_cu  vf_cu, p_cu    
auxOut.Yss        = [Css   Rss    PAIss     Wss   lss    OUTPUTss   MCss   EVFss   VFss   P     ]'; 
auxOut.transformY = [ones(1,9) ones(1,40)];
auxOut.labely     = [{'$c_t$'},{'$r_t$'},{'$\pi_t$'},{'$w_t$'},{'$l_t$'},{'$output_t$'},{'$mc_t$'},{'$evf_t$'},{'$mvf_t$'},nameP];

% Additional output
auxOut.params         = params;
auxOut.params.PHIzero = PHIzero;
auxOut.params.Kss     = Kss;
auxOut.params.OUTPUTss= OUTPUTss;
auxOut.params.Css     = Css; 
auxOut.params.AA      = AA;
auxOut.params.Rss     = Rss;
auxOut.params.MCss    = MCss;
auxOut.params.Wss     = Wss;

%% non-model specific part of output
%The transformed levels of controls
auxOut.yssTrans   = zeros(ny,1);
for i=1:ny
    if auxOut.transformY(i) == 1
        % Log-approximation
        auxOut.yssTrans(i,1) = log(auxOut.Yss(i,1));
    elseif auxOut.transformY(i) == 2
        % Logistic transformation
         auxOut.yssTrans(i,1) = -log(1/auxOut.Yss(i,1)-1);
    elseif auxOut.transformY(i) == 0
        % No transformation    
        auxOut.yssTrans(i,1) = auxOut.Yss(i,1);
    end
end

%The transformed levels of states
auxOut.xssTrans   = zeros(nx,1);
for i=1:nx
    if auxOut.transformX(i) == 1
        % Log-approximation
        auxOut.xssTrans(i,1) = log(auxOut.Xss(i,1));
    elseif auxOut.transformX(i) == 2
        % Logistic transformation
         auxOut.xssTrans(i,1) = -log(1/auxOut.Xss(i,1)-1);
    elseif auxOut.transformX(i) == 0
        % No transformation    
        auxOut.xssTrans(i,1) = auxOut.Xss(i,1);
    end
end
auxOut.ny         = ny;
auxOut.nx         = nx;
auxOut.nx1        = nx1;
auxOut.ne         = ne;
auxOut.mx         = mx;
auxOut.myx        = myx;
auxOut.eta        = eta;
auxOut.momEps     = momEps;




end